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The analysis, understanding, and comparison of network structures is a prominent topic not 
only in computer science theory but also in a diverse set of application-oriented fields. For in¬ 
stance, social media sites use tools from this domain to understand large-scale structures arising 
from localized connections between users, while computer graphics software applies the same 
methodology to understanding the shapes of meshed surfaces. Further applications of this toolkit 
range from designing random walk models using Markov chains [Lov93] to characterizing pat¬ 
terns in molecules, neural networks, and food webs [CCK07], 

The obvious structure for representing a network is a graph, given by a collection of nodes 
connected by edges (optionally associated with distances or weights). The analysis of graphs as 
discrete objects is a classical topic predating the development of computer technology. Algorithms 
for finding shortest paths, routing flow, finding spanning subtrees, identifying cycles, and so on 
are well-understood, although many simple problems about graphs remain unsolved. 

Some more modern branches of graph theory, however, treat graphs themselves as mutable 
objects that should be understood in a continuous or probabilistic context. Graphs from many 
sources are subject to discrete and continuous changes: social networks grow and evolve as users 
"friend" and "un-friend" each other, links can lengthen or strengthen, and so on. Due to this 
mutability as well as the possibility of error in the collection of graph data, multiscale methods 
characterizing graphs at different levels of detail are helpful for understanding broad or changing 
patterns in connectivity. Parameters such as scale and sensitivity, however, can be understood 
as continuously-varying, and thus research carrying out such "semi-discrete" analysis with rigor 
must combine insights from discrete graph theory and real analysis. 

Of special interest are methods making use of partial differential equations (PDEs) and spec¬ 
tral analysis. PDE methodologies model how physical interactions on a local scale can charac¬ 
terize global behavior of waves, heat, and other phenomena on assorted media and domains. A 
fundamental and predictable theme in PDE analysis is that the operators used to construct PDEs 
encode information about the structure of the domains with which they are associated. For exam¬ 
ple, the differential building blocks of the heat equation determine how a given surface conducts 
heat. Less obviously, the De Rham cohomology from differential topology demonstrates that the 
topology of a smooth surface can be understood by analyzing applications of its "first derivative" 
operator d [MT97], 

While first-order PDEs largely can be solved exactly, understanding the second-order Lapla- 
cian operator leads to techniques for analyzing functions on a surface and the geometry of the sur¬ 
face itself [Ros97], Applications in geometry processing discretize these equations using finite ele¬ 
ments, discrete exterior calculus (DEC), and other approximations of the Laplacian for per-vertex 
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signals on meshed surfaces [Hir03, RWSN09]. Many such discretizations exhibit convergence in 
the limit of mesh refinement, providing a connection to the theory of smooth surfaces, and in the 
case of DEC discrete theorems provide exact descriptions of behavior without any limiting pro¬ 
cess. Regardless, with a discrete surface Laplacian in place, tools for matching, segmentation, and 
other tasks can be designed using solutions to the heat and wave equations, among other model 
PDEs [Broil], 

PDE-based analysis of discrete surfaces has the advantage that it can draw intuition from par¬ 
allel constructions in differential geometry. Differential graph analysis, on the other hand, requires 
the development of a new framework to understand diffusion, oscillation, and other phenomena 
as they might occur on graphs, which may not be physically realizable. Instead, we define a Lapla¬ 
cian for functions on a graph with particular manipulations, constructions, and applications in 
mind. Of particular success has been the field of "spectral graph theory," which involves study¬ 
ing eigenvalues of matrix-based graph Laplacians and connecting their properties to those of the 
graph itself [Chu97]. 

Spectral graph theory still is a fairly discrete area of graph analysis in that Laplacians in this 
domain are sparse matrices and their resulting spectra are finite. In contrast, solving PDEs like 
the heat and wave equations on a graph involve continuous derivatives in a time variable. In this 
nascent area of research, there is less consensus on the choice of appropriate "differential" opera¬ 
tors and treatments of the time variable to yield meaningful solutions. Such solutions, however, 
can yield insight into how signals might propagate from node to node of a graph through its 
edges. 

Here, we will discuss and contrast approaches proposed for graph adaptations of three model 
PDEs: the Poisson, wave, and heat equations (introduced in [CYOO, FT04b, CCK07], resp.). In par¬ 
ticular, we will show how the choice of a graph Laplacian operator complements the solutions and 
properties of each model equation. Such a direct comparison motivates discussion of modeling 
more complex PDEs on graphs and potential directions for future research. 

1 Basic PDE Analysis 

Here we introduce the model PDEs we eventually will consider on graphs. We focus on the defini¬ 
tion and structure of each of these equations rather than its derivation from physical principles; we 
refer the reader to any basic PDE textbook, such as [Str08], for a physically-motivated discussion. 

In general, we will consider a partial differential equation in a scalar function u(x) to be any 
equation of the form F({D n u : n E Z + },x) = G(x), where D n u denotes the set of nth derivatives 
of u. To simplify matters, we will consider only the case where F is linear in u and its derivatives, 
thus denoting a linear PDE. We take u to be a sufficiently differentiable function u : O —> R for 
some domain Q; optionally we can include an additional "time" variable t E R + in which case 
u : Q x R + —>■ R. We remain purposefully vague in this section with respect to the choice of Q; 
of course one can think of Q C R !I for the usual discussion of model PDEs, but we eventually will 
take Q to be a graph. 

1.1 The Laplacian 

The "method of characteristics" readily provides solutions to most model first-order PDEs in u, 
but including a second-order term is sufficient to yield diverse and nontrivial behavior. To this 
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end, we introduce a linear Laplacian operator A, which evaluates the total second derivative of u; 
for instance, the Laplacian of u : Q C R” —y R takes the following form: 



d 2 u 
dx 2 


( 1 ) 


Symbolically we can factor the Laplacian operator from the right hand side, and thus we consider 
A as a functional A : C°°(Q) —* C°°(0). In the case of Q C R" above we thus write 



<P_ 

dx 2 


( 2 ) 


In functional analysis, the Laplacian is the canonical example of a compact self-adjoint linear 
operator. Specifically, the operator satisfies the following properties: 

Linear For u, v G C°°(Q) and ci,C 2 € R, A(ci« + civ) = c\Au + C 2 AU. 

Compact A is bounded, in that for appropriate functional norm || • || there exists some M > 0 
satisfying ||Au|| < M||w|| for all u G C°°(0). Furthermore, A applied to a sufficiently small 
open neighborhood of the function u = 0 using the topology induced by the norm || • || 
yields a subset of a compact set in the same functional space. 

Self-adjoint Ifwetakeu,u G C°°(n),then (Aw, v) = (u,Av) for the functional inner product (■, ■). 

Many of these properties are difficult to verify for functional operators on O C R" but become 
trivial for operators such as the graph Laplacian, which operate on finite-dimensional spaces. 
In this case, self-adjoint operators are given by symmetric matrices M G R” xn (thus satisfying 
M = M t ) since we can write (Mv) ■ tv = ( Mv) r w = v T (M T to) = v ■ (M t zh). 

To distinguish the Laplacian from other elliptic operators, we will give it one more property 
that will make it straightforward to generalize proofs of certain analytical theorems: 

Categorizes extrema If x G O is a local minimum of u : O —> R, then [Ah] (jc) < 0. 

Flere, we will assume that the idea of a "local minimum" is adapted to the domain Q in question. 
For instance, if Q C R” and u G C°°(n), then x will be a local minimum when u(x) < u(y) for all 
y satisfying ||x — y|| < e for some e > 0. On a graph, a local minimum of a function on vertices 
will be a vertex whose associated value is less than or equal to the values at its neighbors. 

One principal reason that self-adjoint linear operators are so appealing is that they admit a full 
set of orthogonal eigenvectors spanning L 2 (Q). 1 In particular, there exist orthonormal eigenfunc¬ 
tions (p\, (p 2 , ■ ■ ■ with corresponding eigenvalues 0 < < A 2 < ... such that / G L 2 (0) can be 

written / = X/ a i<pi- 

For such a decomposition to exist, in many cases we must introduce boundary conditions. 
Thus, if desired or necessary we decompose Q into a disjoint boundary and interior as Q = dCl U 
Q. In this case, we will consider only eigenfunctions cpj satisfying A/, = A j(pj in Q and <j>/|an = 0. 
This restriction corresponds to Dirichlet boundary conditions; we will not consider derivative- 
based Neumann boundary conditions or others here, although most of the frameworks we discuss 
will behave equally well in this case. 

1 This discussion is again cavalier with respect to the exact properties of Q for such spectral theory to hold. Since we 

principally consider graphs and compact subsets of R”, we will not delve into the particulars of this theory. 
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1.2 The Elliptic Case: Laplace and Poisson Equations 

The simplest possible PDE involving the Laplacian is the Laplace equation, given by: 


Am = 0 (3) 

This PDE is slightly generalized by the Poisson equation, which prescribes nonzero values for A u: 

A u=f (4) 

The Poisson and Laplace equations are elliptic, since the highest-order derivatives appear in the 
elliptic operator A. 

These equations appear in many situations, from Maxwell's equations for electrostatics and 
the governing equations of Brownian motion to byproducts of the Cauchy-Riemann equations 
defining complex holomorphic functions. In discrete form, they are used in image smoothing 
because they are critical points of the functional f || V7|| 2 dA for image I. 

Note that if we write / = Yi a, (pi for Laplace eigenfunctions (pi and similarly decompose u = 
Yi (pi, then (4) becomes Yi Afifi = L; a i<pu and by orthogonality of the (pi's we have bj = Va,; of 
course if A; = 0 then (pi satisfies the Laplace equation and can be added with a free coefficient. In 
this way, solving the Poisson equation is trivial in the Laplace eigenfunction basis. 

1.3 The Parabolic Case: The Heat Equation 

By introducing a time variable t, we can model PDEs that evolve from some starting condition. 
For instance, consider the time evolution of u given by 

M t = —Am (5) 

where we denote u t = where u satisfies m| t=o = / G C°°(Q) and we constrain u | r)f ) = 0 Vf > 0. 

This PDE is known as the heat equation, because it models the diffusion of / over Q for time t. It 
is a model parabolic equation, because the matrix of high-order derivatives is positive semi-definite 
(only first derivatives appear in t). Similar to the Poisson case, if we write / = Yi a i ( Pi then it is 
easy to see that u = Y a i e ~ X ' t( Pi satisfies our conditions. 

Despite the generality of our choice of O, we actually can prove some facts about solutions of 
the heat equation: 

Lemma 1 (Weak Maximum/Minimum Principle). Suppose u G C°°(Q) satisfies u t > —Am with 
u lan > 0 Vf G [0, T) and u\t=o > 0. Then, u > 0 Vf G [0, T], 

Proof. Our proof follows Theorem 3.1 of [CCK07] but is a standard proof in analysis. Take To G 
(0, T) and £ > 0, and define v = u + ef G C°°(0 x [0, T)). Suppose v has a minimum at (xq, to) G 
Q x [0, To] with xo G Q and fo G (0, To]. Then, we must have 

v t (xo,t 0 ) = u f (x 0 ,f 0 ) +£ < 0 (6) 

because (xo, fo) is a stationary point of v; we only require v t < 0 rather than v t = 0 in case fo = To- 
By the "categorizes extrema" property, we also have 

Av(x 0 , t 0 ) = Am(x 0 , t 0 ) < 0 (7) 
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Combining (6) and (7) yields 


v t {xo,t 0 ) + Av(x 0 ,t 0 ) < o. 


( 8 ) 


This contradicts the relationship 

v t (x 0 ,t 0 ) + A v(x 0 ,t 0 ) = [u t {xo,t 0 ) +Au(x 0 ,t 0 )] +£ > £ > 0 (9) 

by the premises of the Lemma. Thus, we must have that the only minima of v occur when to = 0 
or Xq G 90. Suppose the minimum occurs at (xq, to). For generic (x, t), we have 

u(x, t ) = v(x, t) — et by definition of v 

> v(xq, to) — et by definition of (xo, to) 

= u(x o, to) + efo — ef again by definition of v 

> e(to — t) because u > 0 on 90 and when t = 0 

> — eT since to, t G [0, T] 

Note that T is fixed and £ > 0 is arbitrary, so in reality this inequality shows u(x, t) > 0, as 
needed. □ 

This principle leads us to the uniqueness of solutions to the heat equation. In particular, if u j 
and «2 both satisfy the heat equation with the same boundary conditions, then H h iff = «i — 112 
satisfies the heat equation with zero initial and boundary conditions. By the maximum principle 
applied to H h iff and — M^iff, we must have U h iff = 0 for all f, as needed. In fact, letting t —> 00 
approaches solutions to the Laplace equation Ah = 0, which thus has the unique solution u = 0 for 
our zero boundary conditions. This observation in turn leads to uniqueness of Poisson equation 
solutions A u = g since two solutions Hi and 112 satisfy A (it \ — 112 ) = 0. 

1.4 The Hyperbolic Case: The Wave Equation 

Our final model second-order PDE is the wave equation, given by 

u tt = -An (10) 

with the same boundary conditions as the heat equation. We also prescribe iq = g at t = 0. 
Once again we can compute solutions in closed form as u = YLi a i cos (y/Xit + bi)(p t , where a,, b, are 
chosen to satisfy the given initial conditions; if A, = 0 then we replace cos(- ■ ■) with a t + bg. It 
is hyperbolic because the time derivative has a different sign than the space derivative. We defer 
discussion of the uniqueness of wave equation solutions until we have developed appropriate 
notions of the gradient of a function, which depends more strongly on the domain Q in question. 

2 Geometric Realizations of Graphs 

We turn now to the case when Q is a graph G = (V,E) for vertices V and undirected edges E. We 
assume \ V\, |E| < 00 and that each edge e G E is associated with a length l e > 0. This description, 
however, is a theoretical construct, so we must describe what it means to diffuse heat or propagate 
waves along a graph; different models of diffusion or differential structures will lead to different 
adaptations of the model PDEs from §1. 
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Symbol 

Description 

Notes 

G = (V,E ) 

Graph 

V is the set of vertices; E C V x V is the set of edges 

le 

Edge length 

Associated with edge e E E 

G 

Geometric realization 

Constructed by gluing intervals C R of length l e 
using the topology of G 

C k {G) 

Differentiable functions 

C k on edge interiors; C° at vertices 

TG 

Tangent bundle 

Union of tangent bundles of edge interiors 

X 

Vector field 

Function G\V —t TG 

v/ 

Gradient 

Vector field on G\V constructed using ID calculus 

V calc ' X 

Calculus divergence 

Scalar function on Q\V 

V 

Vertex measure 

Discrete measure on Q with V(u) = 1 Vu E V 

£ 

Edge measure 

Lebesgue edge measure on G) £{e) = Z e Ve E E 

dY 

Integrating factor 

Written a dV + /I d£ for integration using V and £ 

dV x 

Divergence 

Integrating factor with V ca i c on E and discrete 
terms on V 

dC f 

Laplacian 

Integrating factor given by —dV^/f 

A E 

Edge-based Laplacian 

Coefficient of d£ in dCf given by — V ca i c • V/ 

Av 

Vertex-based Laplacian 

Coefficient of dV in dCf given by n ■ V/ 

A 

Normalized adjacency matrix 

The adjacency matrix of G with rows normalized to 
sum to 1 

A h 

Set expansion 

For A C G, A h = {x E G ■ dist(x. A) < h} 


Figure 1: Notation for our discussion of [FT04b], 


We commence our discussion with the most intuitively-appealing construction, which is able 
to apply analytical results directly but in the end suffers from considerable mathematical chal¬ 
lenges. Our development will follow that of [FT04b]; we summarize notation in Figure 1. 

2.1 Preliminaries 

By far the most common visualization of a graph assigns a point in 1R" to each vertex in V and 
connects the pairs in E using line segments or curves of length l e . Note that such an embedding is 
a construction rather than a fundamental property of graphs. For instance, there is no straightfor¬ 
ward reason why a social network graph should be embeddable in a low-dimensional space in a 
meaningful way. Even so, much of our intuition for graph theory comes from such a construction. 

Formalizing this notion, we define the geometric realization of G as a topological object Q con¬ 
sisting of a closed interval of length l e for each edge e E E whose endpoints are identified with 
those of other segments using the connectivity of G. We identify each e E E with its corresponding 
segment within Q and each v E V with the corresponding point; we will use "vertex" and "edge" 
to refer to either type interchangeably when the difference is unambiguous. 

The maximal open interval contained within each edge e of Q is the interior of e. Note that edge 
interiors are simply line segments and as such can be modeled to conduct heat or waves identi¬ 
cally to segments of R. Thus, the primary challenge of this setup is the treatment of differential 
operators near the points in V, which are nonmanifold when vertices have degree 7 ^ 2. Thus, 
to take advantage of this structure, before developing PDE theories for Q , we develop its global 
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differential structure. 

We first define coordinate-free differentiation of scalar functions on Q. Take TQ to be the "tan¬ 
gent bundle" of Q given by the union of the tangent bundles of the edge interiors, isomorphic to 
(5\V)xR. Then, for functions/ G C k (G) = C k (G\V) nC°(£) and vector fields X G C k (G\V,TG) 
we can define the gradient V/ G C k 1 (Q\V,TQ) and calculus divergence V ca i c ■ X G C k 1 (^\V r ) by 
applying the usual definitions of gradient and divergence within edge interiors. 

It is tempting to define the Laplacian of a function on Q in a similar fashion to complete our 
arsenal of differential graph operators. Unfortunately, doing so trivializes the heat and wave equa¬ 
tions, as we have not specified any way for information to be communicated across the vertices. 
Thus, our development must be considerably more subtle. 

2.2 Integration and Divergence 

Our adaptation of differential operators thus far has considered the one-dimensional structures- 
or edges-of Q but has neglected the zero-dimensional vertices. Interestingly, vertices of Q appear 
in our special treatment of the non-manifold points in Q in a form similar to boundary conditions 
of PDEs on intervals in R. This observation makes sense in light of the fact that we can consider 
diffusion and wave propagation along Q to satisfy the usual equations within edge interiors, with 
boundary conditions that their values match at vertices. 

To facilitate subsequent discussion, we introduce two measures on Q\ 

• A discrete vertex measure V satisfying V(v) = 1 Vu G V. 

• A edge measure £ coinciding with the usual Lebesgue measure on edge interiors. 

Note that this definition is somewhat more restrictive than that in [FT04b], but our results gener¬ 
alize easily to their more flexible setting. We denote integration against these two measures using 
f ■ ■ ■ dV and f ■ ■ ■ d£, resp. 

Relating integrals against V and £ will provide a compact methodology for expressing our 
primary results about PDEs on Q. Thus, we introduce integrating factors dT = a dV + f d£, for 
«:V->R and / G Li(G\V)', we may introduce continuity or differentiability conditions on f> as 
need be. Integration of functions on Q is carried out as: 

fdl = J fa dV + f fp d£ (11) 

As an initial example of the types of operations we can handle using this setup, we develop a 
notion of vector field divergence on Q satisfying a theorem analogous to the divergence theorem 
on R". In general, note that our "calculus divergence" is somewhat unsatisfactory in that it yields 
values only on Q\V. While this restriction might be acceptable for directional vector fields, which 
most naturally are associated with points on edges, divergence may be measurable at vertices as 
well as on edges. 

We have two potential strategies for relating a given test function g G C°°(G) and the di¬ 
vergence of a vector field X: integrating the product y V • X or integrating the product —X ■ Vy. 
Applying intuition from integration by parts, these two should be equivalent for the proper choice 
of V. Unfortunately, carrying out the first expression with the calculus divergence does not yield 
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a divergence theorem for graphs, so we explore the second. Define a functional T>x ■ C°°(G) —> 1R 
given by 

V x {g) = -t X-VgdS (12) 

JQ\V 

We prove a useful formula for such an operator: 

Lemma 2 ([FT04b], Proposition 2.16). Suppose we define e to be the unit vector from u to v paralleling 
edge e = (u,v) E E in TQ\ e and take 

C n-X)(u) = £ ?- x \e(v). (13) 

e=(v,u) 6E 

Then, for Xe C^TQ) and g e C°°(£), 

Vx(g) = j locale • X)g d£-(n • X)g dV ], (14) 

where zve assume X can be extended to de = {u, v} by continuity. 

Proof. By Stokes' Theorem, if we restrict to a single edge e = (u, v) we clearly have the relation 


j V calc -Xd£ = e- (X\ e (v) - X \ e («)). 


where the second term represents a (discrete) integral about the boundary of e. Summing over all 
edges in Q yields 

/ V ca i c -Xd£ = f n-XdV 
Jg Jg 

Consider the vector field Xg. Applying this formula yields: 
j Vcalc • (Xg) d£ = J^h ■ (Xg) dV 

= / (h ■ X)g dV by linearity of the dot product 


By the chain rule we can write V ca i c • (Xg) = ( V ca i c ■ X)g + X ■ Vg. Substituting this expression 
into the integral above and reordering terms yields the desired result. □ 


Thus, if we define the divergence of Q to be the integrating factor 

dVx = (Vcalc • x)d£ - (fi ■ X)dV, 

then we satisfy the formula 


gdV } 


X ■ Vg d£ = 0. 


(15) 


(16) 


This formula serves as our divergence theorem for geometric realizations of graphs. In particular, 
if we take g = 1 then we find 

f dV x = 0, (17) 


as expected since the graphs we consider satisfy dQ = 0. [FT04b] shows how to generalize this 
result when a subset of V is marked as dQ. 
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2.3 Geometric Graph Laplacian 

We now can define the Laplacian of a graph's geometric realization: 

dC.f = — dV^f (18) 

Note that in classical notation this definition coincides with the relation A/ = —V • (V/), and that 
our choice of Laplacians parallels that of [FT04b], We define "edge-based" and "vertex-based" 
Laplacians A E and Ay by requiring dCf = A F d£ + AydV; substituting (15) yields: 

A E f = -V calc • V/ (19) 

Ay/ = n ■ V/ (20) 

Note that a straightforward consequence of Lemma 2 is the relationship 

Jgd£ f = J V/ • Vg d£ = j fdC g (21) 

for sufficiently differentiable f,g. It is also easy to see that A E / = 0 implies that / must be 
edgewise-linear; in this case the vertex-based Laplacian Ay coincides with proposed discrete 
Laplacian matrices for graphs. 

Before modeling with dCf, we wish to understand the structure of its eigenspaces. Of course, 
d£ f encodes two separate operators A/- and Ay, whose eigenstructures we must understand sep¬ 
arately. Considering Ay or A E alone is insufficient, since the former operates only at vertices and 
the latter on edge interiors; instead, we will find that appropriate eigenvalue problems couple the 
two operators, using Ay to put boundary conditions on A E . 

Note that Q\V is a collection of open intervals, which clearly does not admit a countable and 
complete set of eigenpairs (/, A ; ). Thus, we must add more conditions to make our problem 
well-posed. We provide an alternative motivation based on a variational problem. 

Suppose we wish to find critical points of the Rayleigh quotient 

I |V/I 2 is 
f I/I 2 i£ 

It is easy to see that uniformly scaling / has no effect on 1Z, so we equivalently find critical points 
of n{f) = f | V/| 2 d£ subject to f |/| 2 d£ = 1. Such a variational problem leads to the Lagrange 
multiplier function: 

A(/,A) = R(/)-a(/ |/| 2 d£-l) = /(I V/| 2 — A|/| 2 ) d£ + A 
Taking the Gateaux derivative with respect to / in the g direction yields: 
dA(f,A; g) = -^A(f + hg, A)| fc=0 

= / ^(\Vf + hX7g\ 2 -A\f + gh\ 2 )\ h= o d£ 

= J(2Vg-\7f-2Afg) d£ 

= 2 f gdC f - 2A J fgd£ by (21) 

= 2 I g(A E / — A/) d£ + 2 j gAy/ dV by definition of dCf 
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So, if / with / |/| 2 d£ is a critical point of 1Z, then it must satisfy both of: 


( 22 ) 

(23) 


A E f = A/ 
Ay/ = 0 


for some A G R. 

This derivation motivates the definition of an edge-based eigenpair of dC f as a pair (/, A) sat¬ 
isfying / G C°°(Q), Ae/ = A/, and Ay/ = 0. Restricting / to an edge interior of Q gives 
f(x) = A cos (cox + B) for parameter x G [0 ,l e \ and co = y/X, since / is an eigenfunction of 
the one-dimensional Laplacian along edges. More globally, our Rayleigh quotient construction al¬ 
lows for a proof similar to that in [EvalO] for PDEs on Euclidean domains of the following spectral 
theorem: 


Lemma 3 ([FT04b], Proposition 3.2). The edge-based Laplacian admits a sequence of eigenpairs (/,, A,) 
satisfying 

• 0 < Ai < A 2 A ■■ ■ 

• ifi} i e z+ forms a complete orthonormal basis for L 2 (Q), where orthonormality is measured with 
respect to dS 

• A; —>■ 00 

In fact, it is straightforward to extend Weyl's Law to show that the number N of eigenvalues 
< A for fixed A grows like y/X. 

Returning to properties specific to the geometric graph Laplacian, we find that the Ay/ = 
0 condition provides a strong characterization of edge-based eigenfunctions. In particular, for 
eigenvalue A > 0, recall that the corresponding eigenfunction / restricted to an edge e = (u, v ) is 
given by f\ e (x) = A cos {cox + B ) for x G [0, l e \, co = \/X, and some A, B G 1R. We evaluate f\ e and 
its derivatives at the endpoints: 


f{u) = f e { 0 ) = AcosB 

f(v ) = f e (Je) = Acos(col e + B) = Acos(col e ) cos(B) — A sin(co! e ) sin(B) 


/'( 0 ) = — AwsinB = —co 


_ , Jip) - cos(o}l e )f(u) 


sm(col e 


by combining the last two expressions 


Summing the final expression around u is exactly the vertex-based Laplacian at u, yielding the 
following condition for values of / at vertex u: 


0 = Ay/|„ = ^ 

e=(u,v)aE 


/(») 


— cos (col e )f(u) 
sin(col e ) 


(24) 


In fact, any co satisfying (24) is an edge-based eigenvalue, since we can use our cosine expression 
for f\ e to fill in the remaining eigenfunction in edge interiors. In other words, we have proved the 
following lemma: 

Lemma 4. (f,co 2 ) with col e kn for any e G E,k G Z is an edge-based eigenpair if and only if co 
satisfies (24) for all u G V and A E f = co 2 f. 
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This expression is encouraging computationally in that it gives a finite system of equations for 
finding edge-based eigenvalues. Unfortunately, the condition is highly nonlinear and nonconvex 
in to, so we cannot in general find the complete spectrum of a graph with arbitrary edge lengths. 

One case in which it is possible to find a closed form for the edge-based spectrum of Q is when 
l e = 1 Ve G E. Then, the denominator of (24) has no dependence on e and can be removed, leaving 
behind a much simpler expression that can be solved for to. 

Lemma 5 ([FT04b], Proposition 3.5). Let A be the adjacency matrix of G after normalizing each row 
to sum to 1. Then, for each eigenvalue A of A, arccos(A) + 2rck and 2n — arccos(A) + 2nk are edge- 
based eigenvalues for k G Z + with the same multiplicities. Furthermore, eigenvalues n + nk occur with 
multiplicity |E| — \ V\; if this value is negative, we subtract the multiplicity from those occurring in A. 

We omit the analytical details of the proof. 

2.4 Geometric Graph Wave Equation 

Having developed differential operators on Q up to second order, we are now equipped to intro¬ 
duce the general wave equation on Q : 

(ctdV + f>d£)u tt = cLZTyVi, (25) 

Here, we make the following assumptions: 

• u G C°(Q x I) with «(■, t) G C 2 (Q) for some interval ICR 

• a, f> G C°(Q) and 7 G C X {Q) with a,f>, 7 > 0 

Note if 7 = 1 then the equation becomes (adV + f>d£)u n = ~dC u , mimicking the classical wave 
equation. Splitting into dV and into d£ terms yields the following system of equations: 


au t t = — 7 A yu (26) 

fu t t = — yA E u + V 7 • Vu (27) 

To prove basic properties of (25), we will examine the energy function given by 

Energy(A;f) = /" ( 7 (Vu) 2 d£ + u 2 (a dV + p d£)) (28) 

J A 

For A C Q define A h = {x G Q : dist(x, A) < h}. The following theorem characterizes the 
evolution of this energy function over time: 

Lemma 6. Suppose 7 < c 2 f> on Q. Then, ifu satisfies (25) then Energy(A ct °;0) > Energy(A-, to). 

The proof of this theorem once again follows that of analogous results from classical PDEs. 
It immediately implies uniqueness of solutions to the wave equation given u and u t at t = 0. It 
also shows that if two solutions agree on some subset of Q, then they are guaranteed to agree 
on a subset of Q shrinking at rate c over time. That is, the wave equation (25) with f > 0 can 
be associated with a speed of propagation c at which information can be communicated along Q; 
this property does not hold for most formulations of the wave equation on graphs. For example. 
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taking /3 = 0 corresponds to an edgewise-linear model of wave propagation, and it can be shown 
that a "hat function" about any vertex in Q can propagate with infinite speed in this case. 

Using arguments similar to our earlier discussion, it is easy to see that when all edges of Q have 
the same length we can find a closed-form solution to the wave equation without considering edge 
interiors. Such a solution technique uses the Chebyshev polynomials in the matrix A to expand 
cos (nA) for n E Z + . 

2.5 Discussion and Applications 

[FT04a, FT04b] provide many theoretical applications of the geometric graph Laplacian and other 
differential structures on Q. Spectral bounds on the diameter of G and related quantities can be 
derived by applying analytical results involving Laplacian eigenvalues. Additional results on 
graph structure can be derived by studying the heat equation with respect to the vertex Laplacian, 
by considering the class of edgewise-linear functions [FT04a]; a more generic heat equation similar 
to (25) is not considered. 

This construction of "calculus on graphs" is very similar to the idea of a "quantum" or "metric" 
graph appearing in the theoretical physics literature paired with Neumann or natural matching 
boundary conditions [GS06]. While some results from either field may be applicable to the other, 
the physical interpretations of PDEs on graphs are different in the two cases. For instance, quan¬ 
tum graphs unsurprisingly involve the wave equations appearing in quantum physics, whereas 
our constructions may be better interpreted as more classical oscillations of a graph realized using 
a series of wires. 

Such geometrizations appear to be useful tools for proving theoretical bounds on graph struc¬ 
tures. In particular, many theorems regarding large-scale behavior of assorted differential equa¬ 
tions carry over directly to the graph case with little to no modification, as we have seen studying 
Laplacian eigenfunctions and the propagation of waves. As a computational tool, however, these 
constructions can be cumbersome. In particular, when edge lengths are nonuniform, we have no 
straightforward characterization of Laplacian eigenfunctions or solutions of the wave equation 
and are forced to resort to approximation by chains of unit edges or finite elements simulation. 

3 Discrete Operators on Graphs 

Rather than posing PDE problems on the continuous topology given by Q , we can attempt to for¬ 
mulate analogous PDE-type equations on the original discrete graph G. In this case, rather than 
applying PDE theory directly, we use it to inspire potential fully-discrete results. These tech¬ 
niques have the advantage that they generally involve finite-dimensional linear algebra in ]Rl y l 
rather than more complex differential theory, even if Laplacians cannot be derived as easily from 
continuous theory. 

3.1 Graph Laplacians 

Taking a reverse perspective for deriving potential graph Laplacians, we can examine patterns 
that appear in discretizations of continuous PDEs, in hopes that the Laplacian matrices from 
these domains have similar properties to the continuous operators that appear in the limit. For 
instance, suppose we approximate a one-dimensional u : R —> ]R with a sequence of values 
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..., U-i, Ho, U \,.... Simple finite differencing approximates du /dx with the sequence differences 
Ui — ii{- 1 - Taking differences of these differences yields a second derivative m,- + i — 2 Uj + iq-_i, or 
convolution with a filter given by ( 1 -2 1 ). A similar approach on an m- dimensional grid 

yields a weight of — n on each central vertex and 1 on each neighbor. This construction suggests a 
"combinatorial" Laplacian for per-vertex functions given as vectors in RlU, given by 

{ deg u if u = v 

— 1 if u and v are adjacent (29) 

0 otherwise 

Similarly, consider random walks along a graph G with transition matrix P. That is P(u,v ) 
is the probability that state u transitions to state v where u and v are adjacent in the graph, and 
otherwise P(u, v) = 0. A positive vector x E RlU summing to 1 can represent the probability that 
a particle is at one of the vertices of G; Px then gives the distribution of potential particle positions 
after one time step. 

One might study stationary distributions of such a Markov process, that is, probability distri¬ 
butions on the vertices of a graph that do not change after a step of the random walk. These 
distributions are given by the equation Px = x, or equivalently (7 — P)x = 0. Notice the form of 
the matrix I — P: 

{ 1 if u = v 

—P(u, v ) if u and v are adjacent (30) 

0 otherwise 

In particular, if a vertex u has equal probability of transitioning to any of its neighbors, this matrix 
equal to T _1 L, where T is the diagonal matrix of vertex degrees. 

Motivated by these two examples and using the notation of [CY00], for a graph G = (V,E) 
with edge weights w e we define the combinatorial Laplacian matrix L of G to be: 

{ d u if u = v 

—W(u,v) ^ 11 an( i v ar e adjacent (31) 

0 otherwise 

with d u = Yj( u , v )eE w (u,v)- Using our above example and analogous appearances of the Laplacian in 
continuous probability theory, we also define a discrete Laplacian A = T~ 1 L, for a diagonal matrix 
T containing the values d u . Neither L nor A is necessarily symmetric, so we define a symmetric 
normalized Laplacian matrix satisfying «£? = T 1/,2 AT~ 1/2 = T~ 1/2 LT~ 1/2 . 

We take S C V to be a subset of vertices of the graph and say a function g E IRU satisfies 
Dirichlet boundary conditions if it is zero on V\S. We denote the Dirichlet eigenvalues Ai < A 2 < 

■ ■ ■ < A s of the graph, for s = |S| to be eigenvalues of 2z?s, the submatrix of S£ corresponding to 
the vertices in S. 

3.2 Green's Functions 

The Green's function of a PDE is its solution operator. Specifically, suppose a PDE on R takes the 
form Lu(x) = f(x) for unknown u and linear differential operator L. Then, the Green's function 
G : R 2 —> R (more precisely a distribution) provides solutions u for any / by convolution: 

u(x) = J G{x,s)f{s) ds (32) 


13 


Integration against G effectively provides an "inverse" for L, although boundary conditions must 
be established to justify such a statement. 

While rigorously establishing that G provides a unique solution to a given PDE can be an ana¬ 
lytical challenge, the construction and manipulations of discrete Green's functions is much simpler. 
In particular, operators like the graph Laplacians from §3.1 are simply matrices in x l y l, so their 
Green's functions intuitively can be computed using matrix inverses. 

One way to understand the construction of discrete Green's functions is through eigen-analysis 
on the matrix «2g. Imitating the proof for geometric Laplacians in §2.3, we can see that the minimal 
eigenvalue Ai of can be computed as: 

. . , ^,i/eSU(5s(/(' 1 ) — f(y)) w (x,y) 

/ Lxesf 2 ( x ) d x 

Here, we define the boundary SS as those vertices in V\S adjacent to vertices in S. If S is connected 
and SS is nonempty, then this expression shows Ai > 0, and thus we can construct an inverse if of 
Jz?S- We consider to be an operator on IRWl, extending to vertices in V\S using zeros. Then, we 
can write a matrix G = T | / 2 6 f 7 | /- providing a discrete Green's function for A; that is, it satisfies 
GA = J on S. 

One can use Green's functions of elliptic operators like the Laplacian to interpret solutions to 
more complex PDEs. For instance, suppose for t > 0 we define the heat kernel = e l ^ s using 
matrix exponentiation. Then, the following heat equation is satisfied: 

j t Mf = -&sf (34) 

The Green's function in this case appears using the relationship if = J 0 °° dt. For this reason, 
discrete Green's functions are of interest as the basic building blocks for studying not only their 
associated operators but also other equations in which they appear. 

While Green's functions may be best described using matrix inverses, algorithmically this de¬ 
scription leads to inefficient techniques: inverses of the sparse matrices Jzfg, L$, and A 5 are not 
guaranteed to be sparse, and we have provided no guarantees on their conditioning. If we ex¬ 
pand the matrix exponential in (34), however, we see that large eigenvalues of .£ are dampened 
quickly as time t progresses, so it is possible to focus computations on eigenvalues of close to 
0 . 

Since is symmetric and positive definite, we can write its eigenfunctions as <pi , <p 2 , ■■■, <^151 
with corresponding eigenvalues 0 < Ai < A 2 < • • • < A 151 . Applying the usual eigenvector matrix 
decomposition element-by-element, we can write 

^(*/y) = ( 35 ) 

Thus, as A —> 00 , the influence of A on £ vanishes. This heuristic argument aside, if we compute 
the whole spectrum of ‘S we can use it to construct exact solutions to vertex-based PDEs on graphs 
explicitly. 

3.3 Solving the Laplace and Poisson Equations on Graphs 

Suppose we wish to find / satisfying A/ = g for some g with nontrivial values on S U SS. We 
can use the eigenvalue constructions of §3.2 to find such a solution explicitly. We first consider the 
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Laplace equation A/ = 0 and then extend our solution to the more general A/ = g: 

Lemma 7 ([CYOO], Theorem 1). Suppose A/ = 0 on S with f(x) = a(x) for x E SS. Then, f satisfies 2 

/( z ) = 4“ 1 / 2 Er^( z ) E d x 1/2 (pi(x)a(x/) (36) 

i i xeS 

(x,y)eE 
y6(5S 

Proof. For convenience, as in [CYOO] we denote / = T 1 / 2 / and note that A/ = 0 if and only if 
Jz?s/ = 0. Thus, it is straightforward to translate between solutions of Laplace's equation with the 
symmetrized Laplacian ££ and those using the discrete Laplacian A. 

Since / satisfies an equation in «£?$, it makes sense to expand it in the eigenfunctions {(pi} as 
/ = £,■ a,(p,, where by orthogonality of the (pi's we have cij = (c pj , /). For convenience, the boundary 
conditions can be encoded by a function fo(x) given by a{x) on SS and 0 otherwise; we will use 
subscripts of S to denote a function restricted to S. 

With the notation in place, it is possible to perform considerable simplification: 

A rn = A i((pi,f) since a, = (c p it f) 

= (J 2 ; ’s<pi,f) since T? s <Pi = A i(p { 

= (^sfi, (/ - /o)) since f 0 = 0 on S 

= {(pi, T~ 1 / 2 LT~ 1/2 (/ — /o)) since = r p- x l 2 \gp- 1 l 2 anc j i s symmetric 
= {(pi, T~ 1/2 L(f - fo)s) by definition of / 

= {T 1/2 (pj, A(/ — fo)s) since A = T^ 1 L and T is symmetric 
= {T 1/2 (p ir — (A/ 0 )s) since A/ = 0 

= - E ( \f d x<pi{x) ■ -J- E (/o( x ) “ /o(y)) I by expanding the inner product 

V Ux (x,y)eE J 

= Y2 E d x 1/2( Pi( x ) cr (y) since /o = 0 in S and fo{x) = a{x) on SS 

.reS (x,y)eE 

yess 


We divide both sides by A, and substitute into the eigenfunction expansion of /: 


/ 


f = E a i(pi = L 


\ 


y. E d x 1/2< Pi( x My ) 

1 ( x,y)eE 

V yess 


<Pi 


Finally, we apply the relationship / = T 1 / 2 /: 

( 

fs{z) = E 


\ 


y E d x 1/2 <pi(*My) 

1 (x,y)eE 
\ yess 


d z 1 / 2 (pi{z) 


2 As an aside, the exponent of is correctly — i rather than | in [CYOO]. 


(37) 


(38) 
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(39) 



E d 

xeS 

(x,y)eE 

!/eSS 


— 1/2 


(*My) 


This final refactorization completes the proof. □ 

Lemma 8 ([CYOO], Theorem 2). Suppose Af = g on S with f\ss = u, where S is a connected subgraph. 
Then, f can be decomposed as a sum f\ + / 2 , where Afi = 0 with f\ss = c and f 2 = Gg for Green's 
function G of A. 

Proof. Take fi= f — fi■ Then, f = fi~\~ fi and Af\ = Af — A /2 = g — g = 0, as desired. □ 

With this characterization of Poisson equation solutions, it is straightforward to write closed- 
form Green's functions and solutions to Laplace and Poisson equations on graphs with regular 
structures [CYOO]. Predictably, Green's functions of grid-structured graphs nearly coincide with 
the usual Discrete Fourier Transform basis. 

3.4 Green's Functions without Boundaries 

Our construction depended on the fact that the set S of vertices on which we computed Green's 
functions satisfied SS 7 ^ 0. When SS = 0, L = L$ no longer is full-rank, so we must deal with this 
case using some specialized constructions. 

It is easy to see that the vector (po given by (po(k) = \/‘Vy, d, comprises the null space of .£. 
Since «£? thus projects (po out of functions on the graph, we simply exclude it when constructing 
the Green's function ( S\ 

= 1 - (p 0 (p<f (40) 

Furthermore we constrain £ /o/ () = 0. Since the matrix (po<po projects onto unit vector (po, these 
constraints simply state that should project out the (po component and act as an inverse to T£ on 
the subspace orthogonal to (po. 

With such a lack of boundary conditions, it is only possible to solve Poisson-type equations up 
to constant shift by <pQ. Parabolic equations such as the heat equation remain largely unchanged 
since they mostly depend on initial conditions. Other problems including many involving Markov 
chains are solved exactly by projecting out (po [AF02]; construction of such problems is out of scope 
of the discussion at hand. 

3.5 Relationship to Geometric Graph Realization 

The development in this section relied fundamentally on the collection of discrete Laplacians de¬ 
fined in §3.1. These Laplacians satisfy the properties stated in §1.1, and thus generic proofs that 
depend only on "Laplacian-like" structures apply equally well to these Laplacians as to the geo¬ 
metric Laplacians introduced in §2.3. 

To align completely with continuous intuition, however, ideally one might expect discrete 
Green's functions and discrete Laplacian eigenfunctions to coincide exactly with their counter¬ 
parts on the geometric graph realization Q restricted to graph vertices. Such optimism likely 
would remain unsatisfied in the most general case: it seems unlikely to be able to find solutions 
to continuous PDEs like the geometric graph wave equation knowing and making use only of 
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solution values at discrete locations. Even so, we can prove that in certain cases the geometric 
and discrete graph Laplacians have more structure in common than the most basic properties of 
Laplacian operators. 

Recall Lemma 5, which characterizes geometric graph eigenvalues in the unit edge length 
case. These geometric eigenvalues were given by {arccos(A) + 27rZ + ,27r — arccos(A) + 27iZ + : 
A E A(A)}, where A represents an eigenvalue of the normalized adjacency matrix A. Writing A 
element-by-element yields 


A(u,v ) 


1 

deg u 

0 


if n and v are adjacent 
otherwise 


(41) 


This matrix has exactly the same structure as the Markov matrix P from §3.1 when the probability 
of a transition from a vertex to any of its neighbors is uniform. Eigenfunctions of P and I — 
P trivially are identical, and thus we can see that in this case eigenfunctions of the discrete 
Laplacian A coincide with geometric Laplacian eigenfunctions restricted to V. Corresponding 
eigenvalues are related using the arccos A nonlinearity, although it is important to remember that 
Q has an infinite spectrum since there are multiple frequencies within a given edge yielding the 
same boundary values. 

Thus, in the equi-length case the discrete Laplacian actually can be used to construct geometric 
PDE solutions. As edge lengths become nonuniform, however, this relationship breaks down, 
mostly due to the nonconstant denominator in (24). 


4 Semi-Discrete PDEs on Graphs 

In §3, we defined discrete Laplacian operators for per-vertex functions on graphs. These operators 
have key properties in common with Laplacians encountered in classical analysis and admit struc¬ 
tures like Green's functions paralleling objects on IR” and on manifolds. Given this relationship, 
it is a reasonable step to consider using discrete Laplacians to model flows on graphs. Intuitively, 
such models yield flows that are spatially discrete along the domain Q but continuous in time t. 

This construction might be considered a "semi-discrete PDE," in the sense that the time vari¬ 
able remains continuous. An alternative viewpoint is that these equations become Rl y l-valued 
systems of ordinary differential equations (ODEs) with a single independent variable t. Both of 
these dual viewpoints are valuable: 

• The PDE standpoint indicates that resulting solutions likely will have structures resembling 
heat flow, wave propagation, and other physical effects modeled using parabolic and hyper¬ 
bolic equations. 

• The mathematics of ODEs provides straightforward characterizations of existence, unique¬ 
ness, and properties of solutions that are stronger than their PDE equivalents. For example, 
the Picard-Lindelof Existence and Uniqueness Thorem for ODEs automatically guarantees 
that discrete graph heat and wave equations have solutions for all time t > 0, effectively 
sidestepping the need for more specialized existence and uniqueness proofs needed for 
many PDEs [HSD12], 

[CCK07] exhaustively lists solution techniques for first- and second-order semi-discrete PDEs 
on graphs including low-order driving terms and other special cases. Here, we will parallel their 
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development in the homogeneous case, characterizing solutions to the model equations of §1. 
Analogous techniques apply to the more general case, although notation and technicalities become 
considerably more complex. 

4.1 oz-PDEs on Graphs 

We will continue to use the same notation as previous sections but briefly connect it to the notation 
used in [CCK07]. There, graph edges are associated with weights to : V x V — > [0, oo); these are 
identical to the weights W{ x ,y) defined earlier, taking W( x ,y) = 0 when (x,y) £ E. The prefix to is 
added to PDE terminology to disambiguate from the continuous case. Vertex degrees d co u are the 
same as our d u , and integration of / : V —> R is given by 

[ fd w =Y, f{u)d u = [ fd u dV 
JG ue y J 

In contrast to [FT04b, FT04a], [CCK07] develops a discrete calculus on graphs without con¬ 
structing a geometric realization. They define a directional derivative [D v f](u) of / : V —> 1R by 
writing differences 

[D*/](k) ee (/OO “/(«)) (42) 

The gradient V/ : V —> IRj y l at u is defined as a vector of D v f values evaluated at u: 

(V/)(«) EE ( D v f(u)) veV (43) 

Second-order PDEs are modeled using the discrete Faplacian A introduced in §3.1. 3 The definition 
of directional differentiation above is consistent with this Faplacian, in the sense that the Faplacian 
looks like summing double applications of D v . This self-consistent system admits analogs of many 
integration-by-parts and related vector calculus identities for discrete graph operators [CB05]. 

4.2 oz-Diffusion Equations 

Although our definition of the Faplacian has changed, the form of the homogeneous m-diffusion 
or heat equation on discrete graphs remains the same: 

Ut = —Am, 

although now we have u e IR |V| a signal only on vertices rather than Q. Expanding the Faplacian 
row-by-row demonstrates that this equation models flows between nodes where conductivity is 
measured using w and the instantaneous flow rate is determined by the difference in value of u 
between adjacent nodes; such models appear in discretized PDEs as well as in models of electrical 
networks. Notice that the proof of Femma 1 applies for this heat equation, implying not only 
existence and uniqueness but also a characterization of the continuum of values appearing in u. 

Employing the notation of §3.2, we will assume that a subset S C Vis specified on which we 
desire to solve the heat equation and that boundary values are prescribed on SS. To do so, we once 
again consider the normalized Faplacian J 2 ? = T 1,/2 AT 'Z 2 restricted to S, yielding operator «5?s 

3 The Laplacian used here differs from that in [CCK07] by a sign. 
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with eigenvalues 0 < Ai < \i < ■ ■ ■ < A|g| and eigenfunctions (pi,(pi,- ■ ■,f\s\■ Since is positive 
definite, we can assume that the (pj's are orthogonal and unit length. 

We can use intuition from the construction of (35) and the fact that e~^ it (pi satisfies the heat 
equation Vz to construct a heat kernel K$: 


|S| 


K s (u,v,t ) = Yj e Ait <Pi( u )<Pi(. v )\ -y- 


(44) 


i=l 


This kernel provides a succinct description of solutions to the graph heat equation: 


Lemma 9 ([CCK07], Theorem 3.4). Suppose SS f 0. Take functions cr : SS x [0, T) —> ]R and f : 
S — y ]R such that a is continuous and L 1 in time t for each vertex. Then, the solution of the heat equation 
u t = — A u with u | t=o = f and u\$s = cr Vf G [0, T) is given by 


u(x,t ) = (Ks(x,-,t),f(-)) s + J (K s (x,-,t- r), B a (-, r)) dr 

where inner products are over graph vertices summed over and 

a(z,t)w {lJ/Z) 


t) — r. 

zeSS 

(y,z)eE 


(45) 


(46) 


Proof Suppose u is such a solution. By completeness of the basis {(pfi, we can choose functions 
aft) = ( T 1/2 u,(pi ) satisfying 

T l /2u \s = J2 a i(t)(pii x ) (47) 

i 

since each (pj satisfies the zero boundary condition. 

Since «£? = T x ^ 2 AT~ X ^ 2 , we can write S^T 1 ^ 2 = T 1 ^ 2 A. Thus, in the style of the proof of 
Lemma 7, we expand the inner product for a,: 

A;fl;(f) = (T 1/2 u(-, t),\j(pj)s by definition of a, 

= (T 1 / 2 zz(-, t),J£s<pi)s since (pi is an eigenfunction of 
= (T 1 / 2 U f,t),J?(p i ) § -(T 1 / 2 U (;t),J?(p i ) S s where we define S = S U 
= (. £’T 1/2 u(-,t),(pi)s ~ (T 1/2 u(-, t),Jf(pj)ss by symmetry of 
= (T 1/2 Auf,t),(pi)g — (T 1/2 u(-,t),J(f(p i )ss since = T 1/2 AT~ 1/2 
= -aft) - (T 1/2 uf,t),^(pi) ss since u satisfies the heat equation 
= -aft) - ( T 1/2 cr(- r t),3?(pi)ss by boundary conditions on u 


If we expand the definition of A row-by-row, it is easy to see that 

a /(*)= E 

(x,y)6E 
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This implies the following form for applications of Jzf: 

JSf/(*) = [T 1/2 AT~ 1/2 ]f(x) by definition of 2z? 


= 74 E 

(x,y)eE 


fix) fill ) \ w (xy) 

J ' ———■ since we must apply powers of T 


= E 


/(*) f(y) 


(x,y)eE \ d * V^V, 


w 


(*-y) 


In particular, taking / = cp, and evaluating on the boundary x E SS yields 

^47) = - E — 7^4^ 

(x,y)e£ V fl * fl y 

since (pj satisfies the zero Dirichlet boundary conditions. 

Thus, returning to our earlier chain of equalities for shows 

, m , U x , v- cr ( Z ' f )4(l/)^(y,z) 

A Mt) = —« f (0 + E -7T— 

ze<5S v4 

(y,z)eE 

This is an inhomogeneous first-order ordinary differential equation in t whose solution is: 

4(y)^(y,z) 

77 


(48) 


(49) 


flf(f) = Ci e~ Ait + e~ Ait E ( v(z,r)e^ T dT 


(50) 


zgiSS 
(y,z)6E 

for constants {c ; }j =1 | S | C 1R. Note fl, (0) = c, in this characterization. 

To choose the constants c u we return to the definition of fl, (0) and evaluate at t = 0: 

a = (T 1/2 u(-, 0),(pi) = (T 1/2 f,(pi)s 

Finally, we substitute this expression into (47) to show: 

l 

“(*/0 = -73=E fl i(0^W 

v i 

i 


(51) 


x/4 / 


E 


c,e A,f + e x,t E ( / cr(z,r)e AiT dT 

zeS S 
(y-z)e E 


4(y)^(y,z) 

\fdy 


(pi(x) by (50) 


E 


\fdx j 


{TWf'Mse-U + e-U E ( <r(z,T)^ dx 


= E x s(^y^)/(y) 

yes 


ze^S 

(y, z)eE 

v- r, r , ^T>(y,z) 

E K s (x,y,t- t) -n- 

z 6 « 4 

dy,z) 6 E 


<Pi(y)W(y, z) 


4>i(x) by (51) 


\fd~y 

dr by definition of 


= y\K(x,y,t)f(y) + / (K s (x / -,t-r),B cr (- / T)) dr, as desired, 

yes 4 
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This argument shows that if a solution exists to the heat equation then it must take this form. 
It is easy to check that this expression does indeed satisfy the heat equation; in fact, this would 
comprise a short but somewhat obtuse proof of the lemma. Nevertheless, it thus represents the 
unique solution to the heat equation with given boundary and initial conditions. □ 

The proof of Lemma 9 is fairly technical, but from a high level it demonstrates the dual 
PDE/ODE approach to understanding flows on graphs. ODE theorems are used to guarantee ex¬ 
istence and simplify convergence since sums are finite. The construction of the solution, however, 
relies on a separation of variables technique reminiscent of that used to solve the heat equation on 
subsets of R". 


4.3 cu-Elastic Equations 

As with the heat equation, the homogeneous m-elastic equation on discrete graphs takes the usual 
form 

u t t = —Au, 

for signals u E III v . Note that this equation when viewed as an ODE models the behavior of a 
network of linked springs connecting vertices on a graph. 

It is straightforward if verbose to parallel the proof of Lemma 9 to construct closed-form so¬ 
lutions to the cc-elastic equation with given boundary and initial conditions. To do so, we will 
construct an additional Dirichlet kernel for this contrasting problem: 


Ws(x,y,t) 


tyo(*)0o(y) +E 




sin 


(\/A 



(52) 


Then, solutions are provided by the following lemma: 

Lemma 10 ([CCK07], Theorem 4.3). Suppose J>S ^ 0. Take functions a : <5S x [0,T) —* R and 
f,g:S—> R such that a is continuous and L 1 in time t for each vertex. Then, the solution of the elastic 
equation u t = — Au with u\t=o = f, Mf|t=o = g, and u\ss = v\/t E [0, T) is given by 


u(x,t ) = {W s (x,-,t),g) s + 





■,t- r), B a (-, r))s dr 


(53) 


Proof sketch. Define aft) analogously to (47). Paralleling the expansion of A,fl,(i) in Lemma 9 
yields a second -order ODE in aft ) as a function of t, with no first-order term. Applying the usual 
inhomogeneous second-order ODE solution formula and substituting the given boundary condi¬ 
tions yields the desired result. □ 


4.4 Analysis and Properties 

The generic proof of the weak maximum principle in Lemma 1 demonstrates how straightforward 
it can be to extend standard analytical proofs to the semi-discrete case. Such structure preservation 
is a common theme in the treatment of discrete PDEs on graphs. We provide two additional 
examples of continuous theorems that naturally extend to the semi-discrete case, one for the heat 
equation and one for the wave equation. 
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Lemma 11 ([CCK07], "Huygens Property" Theorem 3.6). Suppose u satisfies the heat equation u t = 
—Au on a graph with S = V. Then, for every t,5 > 0 we have 

u(x,t + 5) = (K(x,-,5),u{-,t)) (54) 

Proof. It is easy to check that the solution u in this case satisfies 


«(*/0 = T=E C «' g Xit( Pi( X ) 

V dx 7 


(55) 


for constants {c;}j =1 jyi. Thus, 


\fd~x i 


u(x, t + 5) = -^=Yj c ie Ai{t+S) (pi(x)b y(55) 


\fd~x i 


^e-^ae-^cpfx) 


y=Ee XiS cpi{x) /^fj(pi(y),-j=Y2 c je by orthogonality of the cpfs 


E e x ' S (pi{x){Jd XJ (pi(y),u(-,t)) by (55) again 


vr x j 

= (K(x,-,S),u(-,t)), as desired 


□ 


Our next property provides a discrete analog of (28), which is also conserved over time. 

Lemma 12 ([CCK07], Theorem 4.6). Suppose 5S 7 ^ 0 and u satisfies the co-elastic equation with u\ss = 
0. Then, the following energy function remains constant with respect to t: 



\ E _( DyU(X,t )) 2 

(x,y)eE,yeS 


d x dVx 


Proof. We show that the derivative of E with respect to t is zero: 


(56) 


£, « - L 


2ufUff + E(Dy« ' DyU t ) 

y 


d x dV x , differentiating under the (discrete) integral 

= j [1 lutUtt + Vu ■ Vw t ] d x dV x by definition of the discrete gradient 
= 2 J u t (u tt + A u)d x dV by discrete integration by parts [CB05] 

= 2 j u t (u tt + A u)d x dV since we specified u\ss = 0 
= 0 since u satisfies the wave equation u tt = —Au in S 
Thus, E does not change as t progresses. □ 


Note that this lemma implies that solutions of the discrete wave equation with prescribed 
boundary values are unique. This is because the difference U\ — 112 satisfies the conditions of 
Lemma 12 with E(0) = 0. Thus, ,] l<)i{u\ — uf) = 0, and they agree at t = 0. 
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5 Conclusion 


The abstract definition of a graph as a collection of vertices and edges admits straightforward 
visual and physical intuition but by design does not have a single application or construction in 
mind. Such generality is exactly what makes graphs versatile data structures, but at the same time 
provides limited guidance for the best adaptations of potentially valuable differential techniques. 
Thus, there remains little consensus regarding the most effective methods for modeling flows, 
waves, and other phenomena as they might propagate along a graph. 

Key to the development of a PDE theory for functions on graphs is the choice of topological 
domains and their accompanying operators. Arguably closest to the classical case, [FT04b] argues 
that the geometric realization Q is the most natural domain. Geometrizations have the advantage 
that they can be realized physically and have differential structure in edge interiors, allowing 
for the application of classical one-dimensional results with boundary conditions coupling vertex 
values. Unfortunately, the transition from the manifold structure of subsets of 1R 1 to the "varifold" 
structure of a graph presents considerable theoretical challenges: Even the choice of Laplacian 
operators is unclear and may require combining two operators like Ag and Ay using integrating 
factors. While abstract theorems characterizing graph diameters and other large scale properties 
can be proven in some cases, much remains unknown about Q and whether it can be used to 
construct practical algorithms or analytical tools. 

More practically implementable is the collection of techniques involving discrete differential 
operators on functions in Rl y l. These operators typically are expressible using sparse matrices that 
easily can be computed and manipulated. In particular, depending on application there are at least 
three discrete Laplacian operators that might be useful, notated here as L, A, and ££. Even more 
operators with similar structure can arise when dealing with Laplacians from geometry or simu¬ 
lation applications such as [Hir03]. Thankfully, although the resulting flows can be different nu¬ 
merically, their broad properties remain the same and even can be simplified from the continuous 
case. For instance, discrete Green's functions are constructed using matrix inversion, sidestepping 
the need for specialized convolution constructions and arguments about convergence [CYOO]. 

Partway between the fully continuous and fully discrete approaches are those methods al¬ 
lowing for semi-discrete treatment of PDEs with continuity in time but discrete values along 
graphs [CCK07], These techniques take advantage of the fact that per-vertex functions on graphs 
admit some amount of "discrete differential" structure relating values at neighboring vertices, 
while allowing continuity in variables that are not naturally discretized by the graph itself such 
as time t. These methods provide a potential compromise in which we can understand flows with 
continuous time, although computationally timestepping and eigenfunction computations will 
make related analytical techniques approximate. 

Regardless of the domain and corresponding operator, many common themes arise when con¬ 
sidering flows on graphs; most of these themes are direct adaptations of concepts from classical 
continuous theory. Most importantly, all three of [FT04b, CYOO, CCK07] make use of operator 
eigenfunctions to construct closed-form solutions to model PDEs. In the case of [FT04b], these 
eigenfunctions exist in theory but are countably infinite, whereas in [CYOO, CCK07] they are finite. 
Regardless, the theory of elliptic or positive definite operators guarantees that they fully charac¬ 
terize the behavior of the associated PDEs, and the straightforward form of the heat and wave 
equations simplifies their usage considerably. Eigenvector computation is a well-studied tech¬ 
nique in numerical analysis, so such closed-form solutions can be used directly in applications 
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studying graphs with discrete operators. 

A less technical common theme is that many of the qualitative properties of model PDEs are 
carried over from one domain to the other with little to no adaptation. Many theorems characteriz¬ 
ing solutions to differential equations can be expressed with little more than the list of properties 
of Laplacians in §1.1. This observation confirms the intuition that the heat and wave equations 
have predictable properties on graphs even if they are not directly realizable in Euclidean space. 

We have focused here on theoretical treatments of PDEs on graphs. Work in this domain is not 
necessarily focused on providing algorithmic tools but can be applied directly and indirectly to 
formulate methods for analyzing networks. For instance, having verified that the discrete graph 
Laplacian encodes aspects of structure in a similar way to smooth Laplacians indicates that geo¬ 
metric methods like [SOG09] can be used for graph matching. More generally, the propagation 
of heat and waves from vertex to vertex can help characterize graph topology at multiple scales 
as t —> oo through the use of straightforward numerical routines rather than specialized discrete 
algorithms. Note that drawbacks similar to [GW92] still apply, in that heat and wave kernel meth¬ 
ods have no hope of distinguishing isospectral graphs. 

Considerable work remains to be done toward understanding the advantages and drawbacks 
of network-based models of flows and differential equations. Even so, the promising initial work 
of [FT04b, CYOO, CCK07] indicates that achieving deep results about these models may be as 
simple as adapting pre-existing proofs and techniques from the continuous domain. 
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